# *********************************************************************** #
#              ** Statistical Methods for Data Science I **               #
#              ** Topic 13 - Parametric Inference - MoG **                #
# *********************************************************************** #

# Homemade FWD sampling ---------------------------------------------------

M   <- 100
# M   <- 1000
out <- rep(NA, M)
for (m in 1:M){
  lab = sample(c("Female", "Male"), size = 1, prob = c(.4,.6))
  if (lab == "Female") 
    out[m] <- rnorm(n = 1, mean = 162, sd = sqrt(5))
  else 
    out[m] <- rnorm(n = 1, mean = 176, sd = sqrt(8))
}
hist(out, main = "", xlab = "")

# Bart’s density ----------------------------------------------------------

########################################################################
# X ~ Bart Simpson's density (a.k.a. as "The Claw")                                
#     f(x) = 0.5*dnorm(0,1) + sum(mu in 0:4){0.1*dnorm((mu/2)-1, 0.1)}
# 
# Mixture of 6 Normal densities with:
# pi = c(0.5,  0.1,  0.1,  0.1,  0.1,  0.1)
# mu = c(0.0, -1.0, -0.5,  0.0,  0.5,  1.0)
# s  = c(1.0,  0.1,  0.1,  0.1,  0.1,  0.1)
########################################################################

# Some colors, just in case...
?rgb
colos <- c(rgb(32/255, 74/255, 135/255, 0.7),
           rgb(204/255, 0, 0, 0.7),
           rgb(200/255, 141/255, 0, 0.7),
           rgb(78/255, 153/255, 6/255, 0.7),
           rgb(32/255, 74/255, 135/255, 0.3),
           rgb(204/255, 0, 0, 0.3),
           rgb(200/255, 141/255, 0, 0.3),
           rgb(78/255, 153/255, 6/255, 0.3))

# Package and function
suppressMessages(require(mixtools, quietly = T))

?rnormmix

n <- 250
XX <- rnormmix(n, 
           lambda = c(0.5, rep(0.1,5)), 
           mu     = c(0, ((0:4)/2)-1), 
           sigma  = c(1, rep(0.1,5)) )

# Make an histogram of the data
hist(XX, prob = T, col = gray(.8), border = NA, xlab = "x",
     main = paste("Data from Bart's density",sep=""),
     sub = paste("n = ", n, sep = ""),
     breaks = 50)
# Show the data points
rug(XX, col = rgb(0,0,0,.5))

# Plot the true density
true.den = function(x) 0.5*dnorm(x, 0, 1) + 
                                0.1*dnorm(x,-1.0, 0.1) + 0.1*dnorm(x, -0.5, 0.1) +
                                0.1*dnorm(x, 0.0, 0.1) + 0.1*dnorm(x,  0.5, 0.1) +
                                0.1*dnorm(x, 1.0, 0.1)
curve(true.den, col = rgb(1,0,0,0.4), lwd = 3, n = 500, add = TRUE)

# Kernel density estimate
lines(density(XX),            col = colos[3], lwd = 3)   # Oversmoothing
lines(density(XX, bw = .08),  col = colos[4], lwd = 3)   # Just Right
lines(density(XX, bw = .008), col = colos[5], lwd = 3)   # Undersmoothing

# Add a legend
legend("topright", c("True","Over", "Just right", "Under"), lwd = 5,
       col = c(rgb(1,0,0,0.4), colos[3], colos[4],colos[5]), cex = 0.8, bty = "n")

# Handmade EM4MoG ---------------------------------------------------------

handmade.em <- function(y, p, mu, sigma, n_iter, plot_flag = T)
{
  # Init / 2 components only
  cols     <- c(rgb(1,0,0,.3), rgb(0,1,0,.3))
  like     <- p[1]*dnorm(y, mu[1], sigma[1]) + p[2]*dnorm(y, mu[2], sigma[2])
  deviance <- -2*sum(log(like))
  res      <- matrix(NA,n_iter + 1, 8)
  res[1,]  <- c(0, p, mu, sigma, deviance)
  for (iter in 1:n_iter) {
    # E step
    d1 <- p[1]*dnorm(y, mu[1], sigma[1])
    d2 <- p[2]*dnorm(y, mu[2], sigma[2])
    r1 <- d1/(d1 + d2)
    r2 <- 1 - r1
    
    # M step
    p[1]     <- mean(r1)
    mu[1]    <- sum(r1*y)/sum(r1)
    sigma[1] <-sqrt( sum(r1*(y^2))/sum(r1) - (mu[1])^2 )
    p[2]     <- 1 - p[1]
    mu[2]    <- sum((r2)*y)/sum((r2))
    sigma[2] <- sqrt(sum(r2*(y^2))/sum(r2) - (mu[2])^2)
    
    # -2 x log-likelihood (a.k.a. deviance)
    like     <- p[1]*dnorm(y, mu[1], sigma[1]) + p[2]*dnorm(y, mu[2], sigma[2])
    deviance <- -2*sum( log(like) )
    
    # Save
    res[iter+1,] <- c(iter, p, mu, sigma, deviance)
    
    # Plot
    if (plot_flag){
      hist(y, prob = T, breaks = 30, col = gray(.8), border = NA, 
           main = "", xlab = paste("EM Iteration: ", iter, "/", n_iter, sep = ""))
      set.seed(123)
      points(jitter(y), rep(0,length(y)), 
             pch = 19, cex = .6, 
             col = cols[ (dnorm(y,mu[1],sigma[1]) > dnorm(y,mu[2],sigma[2])) + 1])
      curve(p[1]*dnorm(x,mu[1],sigma[1]) + p[2]*dnorm(x,mu[2],sigma[2]),
            lwd = 4, col = rgb(0,0,0,.5), add = TRUE)
      Sys.sleep(1.5)
    }
  }
  res <- data.frame(res)
  names(res) <- c("iteration","p1","p2","mu1","mu2","sigma1","sigma2","deviance")
  out <- list(parameters = c(p = p, mu = mu, sigma = sigma), deviance = deviance, res = res)
  return(out)
}

data("faithful")
?faithful
hem_fit <- handmade.em(faithful$waiting, 
                       p      = c(.5,.5), 
                       mu     = c(45,55), 
                       sigma  = c(8,8), 
                       n_iter = 20)

round( hem_fit$parameters, 3 )

##     p1     p2    mu1    mu2 sigma1 sigma2 
##  0.360  0.640 54.596 80.079  5.855  5.880

hem_fit$deviance

## [1] 2068.005

# Handmade EM: Different Initial Values -----------------------------------

require(ggplot2)

## Loading required package: ggplot2

require(gridExtra)

## Loading required package: gridExtra

plotConvMC <- function(df, title = NULL)
{
  G      <- (ncol(df)-2)/3
  df$rep <- as.factor(df$rep)
  graf   <- vector("list", ncol(df)-2)
  for (j in (2:(ncol(df)-1))) {
    grafj <- ggplot(df) + geom_line(aes_string(df[,1],df[,j], 
                                               color = df[,ncol(df)])) +
      xlab("iteration") + ylab(names(df[j])) + theme(legend.position = "none")
    graf[[j-1]] <- grafj
  }
  do.call("grid.arrange", c(graf, ncol = 3, top = title))
}


D.em <- NULL
set.seed(1234)
for (m in (1:10)) {
  p1 <- runif(1,0.1,0.9)
  df.em <- handmade.em(faithful$waiting, 
                       p  = c(p1, 1 - p1), 
                       mu = rnorm(2, 70, 15), sigma = rlnorm(2, 2, 0.7),
                       n_iter = 50, plot_flag = F)$res
  df.em$rep <- m
  D.em <- rbind(D.em,df.em)
}
plotConvMC(D.em)

# Up to some permutation (labels switching), all the runs converge to the same solution.
# A very poor initial guess may lead to a very poor convergence of EM.

handmade.em(faithful$waiting, p = c(0.5, 0.5), mu = c(30, 40), sigma = c(2, 10), n_iter = 5)

## $parameters
##     p1     p2    mu1    mu2 sigma1 sigma2 
##    NaN    NaN    NaN    NaN    NaN    NaN 
## 
## $deviance
## [1] NaN
## 
## $res
##   iteration           p1        p2      mu1      mu2       sigma1   sigma2
## 1         0 5.000000e-01 0.5000000 30.00000 40.00000 2.000000e+00 10.00000
## 2         1 1.290643e-11 1.0000000 43.00624 70.89706 1.130252e-01 13.56996
## 3         2 4.706595e-11 1.0000000 43.00000 70.89706 8.259062e-07 13.56996
## 4         3 2.337458e-05 0.9999766 43.00000 70.89771 8.259062e-07 13.56945
## 5         4 3.675314e-03 0.9963247 43.00000 70.99997 0.000000e+00 13.48858
## 6         5          NaN       NaN      NaN      NaN          NaN      NaN
##   deviance
## 1 5227.041
## 2 2190.578
## 3 2190.565
## 4 2174.461
## 5     -Inf
## 6      NaN

# In this example, one of the variances converges essentially to 0...

# mixtools ----------------------------------------------------------------

?normalmixEM
mt_fit = normalmixEM(faithful$waiting)

## number of iterations= 24

summary(mt_fit)

## summary of normalmixEM object:
##           comp 1    comp 2
## lambda  0.360887  0.639113
## mu     54.614891 80.091091
## sigma   5.871244  5.867717
## loglik at estimate:  -1034.002

names(mt_fit)

## [1] "x"          "lambda"     "mu"         "sigma"      "loglik"    
## [6] "posterior"  "all.loglik" "restarts"   "ft"

plot(mt_fit)

class(mt_fit)

## [1] "mixEM"

?plot.mixEM
plot(mt_fit, density = TRUE, w = 1.1)

# Classes
(mt_fit$posterior[,1] < 0.5) + 1

##   [1] 2 1 2 1 2 1 2 2 1 2 1 2 2 1 2 1 1 2 1 2 1 1 2 2 2 2 1 2 2 2 2 2 1 2 2 1 1
##  [38] 2 1 2 2 1 2 1 2 2 1 1 2 1 2 2 1 2 1 2 2 1 2 2 1 2 1 2 1 2 2 2 1 2 2 1 2 2
##  [75] 1 2 1 2 2 2 2 2 2 1 2 2 2 2 1 2 1 2 1 2 1 2 2 2 1 2 1 2 1 2 2 1 2 1 2 2 2
## [112] 1 2 2 1 2 1 2 1 2 1 2 2 1 2 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 1 2 2 2 1 2 1
## [149] 2 1 2 2 1 2 2 2 2 2 1 2 1 2 1 2 1 2 1 2 1 2 1 1 2 2 2 2 2 1 2 2 1 2 2 2 1
## [186] 2 2 1 2 1 2 1 2 2 2 2 2 2 1 2 1 2 2 1 2 1 2 2 1 2 2 2 1 2 1 2 1 2 1 2 1 2
## [223] 1 2 2 2 2 2 2 2 2 1 2 1 2 1 1 2 2 1 2 1 2 1 2 2 1 2 2 2 1 2 2 2 2 2 2 2 1
## [260] 2 2 2 1 2 1 1 2 2 1 2 1 2

# To be compared with its "extreme" version: k-means (with 2 centers)
?kmeans
km_fit <- kmeans( faithful$waiting, centers = c(45, 55) )
names(km_fit)

## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

# Compute the proportion, center and standard deviation for each cluster,
list(p     = km_fit$size/sum(km_fit$size), 
     mu    = as.vector(km_fit$centers), 
     sigma = sqrt(km_fit$withinss/km_fit$size)
     )

## $p
## [1] 0.3676471 0.6323529
## 
## $mu
## [1] 54.75000 80.28488
## 
## $sigma
## [1] 5.865791 5.610953

#...compare...
list(p = mt_fit$lambda, mu = mt_fit$mu, sigma = mt_fit$sigma)

## $p
## [1] 0.3608869 0.6391131
## 
## $mu
## [1] 54.61489 80.09109
## 
## $sigma
## [1] 5.871244 5.867717

# Model Based :: MoG ------------------------------------------------------

## Geyser Data ---
suppressMessages(require(plotly, quietly = T))
?MASS::geyser
dd <- MASS::geyser
kd <- with(MASS::geyser, MASS::kde2d(duration, waiting, n = 50))
p <- plot_ly(x = kd$x, y = kd$y, z = kd$z) %>% add_surface()
p

# Mixture of Gaussian
suppressMessages(require(mclust, quietly = T))
?mclust

mod1 = Mclust(dd)
mod1

## 'Mclust' model object: (VVI,4) 
## 
## Available components: 
##  [1] "call"           "data"           "modelName"      "n"             
##  [5] "d"              "G"              "BIC"            "loglik"        
##  [9] "df"             "bic"            "icl"            "hypvol"        
## [13] "parameters"     "z"              "classification" "uncertainty"

summary(mod1)

## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVI (diagonal, varying volume and shape) model with 4 components: 
## 
##  log-likelihood   n df       BIC       ICL
##        -1330.13 299 19 -2768.568 -2798.746
## 
## Clustering table:
##  1  2  3  4 
## 90 17 98 94

mod1$classification

##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   1   4   3   1   1   4   3   1   4   3   4   3   4   3   1   4   3   4   3   1 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
##   4   3   4   3   4   3   2   1   1   1   1   1   4   3   2   3   4   3   4   3 
##  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
##   4   3   4   3   4   3   4   3   4   3   2   3   4   3   4   3   1   2   1   1 
##  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
##   4   3   4   3   4   3   4   3   1   4   3   4   3   1   1   4   3   1   1   1 
##  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
##   1   4   3   2   1   4   3   4   3   4   3   4   3   4   3   4   3   4   3   4 
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 
##   3   4   3   4   3   4   3   4   3   2   3   1   4   3   4   3   4   3   4   3 
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 
##   4   3   1   1   4   3   1   1   1   1   1   1   2   3   1   1   1   1   4   3 
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 
##   1   1   1   1   1   1   4   3   2   3   4   3   4   3   4   3   4   1   1   1 
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 
##   1   1   1   4   3   4   3   4   3   2   3   2   1   4   3   2   3   4   3   1 
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 
##   4   3   4   3   1   1   1   4   3   4   3   4   3   4   3   2   1   4   3   4 
## 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 
##   3   4   3   1   4   1   1   2   1   1   1   4   3   4   3   4   3   2   1   1 
## 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 
##   4   3   1   1   1   1   1   1   4   3   4   3   4   3   1   1   1   4   3   1 
## 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 
##   4   3   2   1   4   3   2   4   3   4   3   1   1   4   3   4   3   1   1   1 
## 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 
##   1   4   3   1   1   4   3   4   3   2   3   1   4   3   4   3   1   1   1   1 
## 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 
##   1   1   1   4   1   4   3   4   3   4   3   4   3   4   3   4   3   1   4

# Default plot
class(mod1)

## [1] "Mclust"

?plot.Mclust

par(mfrow = c(2,2))
plot(mod1, what = "uncertainty")
plot(mod1, what = "classification")
plot(mod1, what = "density")
plot(mod1, what = "BIC")

par(mfrow = c(1,1))

## Circle ---

# Noiseless
set.seed(123)
n      <- 500
radius <- 2
angle  <- runif(n)*pi*2
x      <- cos(angle)*radius
y      <- sin(angle)*radius

# Jitter around
SD <- .15
xn <- x + rnorm(n, sd = SD)
yn <- y + rnorm(n, sd = SD)
XX <- cbind(xn, yn)
plot(x,y, pch = 19, col = "black", asp = 1, cex = .5)
points(xn, yn, pch = 19, col = "red", cex = .5)

mod2 = Mclust(cbind(xn,yn))
mod2

## 'Mclust' model object: (EEV,8) 
## 
## Available components: 
##  [1] "call"           "data"           "modelName"      "n"             
##  [5] "d"              "G"              "BIC"            "loglik"        
##  [9] "df"             "bic"            "icl"            "hypvol"        
## [13] "parameters"     "z"              "classification" "uncertainty"

summary(mod2)

## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust EEV (ellipsoidal, equal volume and shape) model with 8 components: 
## 
##  log-likelihood   n df       BIC       ICL
##        -1013.33 500 33 -2231.742 -2373.283
## 
## Clustering table:
##  1  2  3  4  5  6  7  8 
## 76 45 76 58 45 69 66 65

par(mfrow = c(2,2))
plot(mod2, what = "uncertainty")
plot(mod2, what = "classification")
plot(mod2, what = "density")
plot(mod2, what = "BIC")

par(mfrow = c(1,1))

# kmeans++ ----------------------------------------------------------------

# k-means++ is a routine to suggest center points before the classical k-means 
# is called. The following lines of code will do that, where X is a matrix of 
# data points, as requested for kmeans, and k the number of centers:

suppressWarnings(require(pracma, quietly = T))
?distmat
?kmeans

kmpp <- function(X, k) {
  n <- nrow(X)
  C <- numeric(k)
  C[1] <- sample(1:n, 1)
  
  for (i in 2:k) {
    dm <- distmat(X, X[C, ])
    pr <- apply(dm, 1, min); pr[C] <- 0
    C[i] <- sample(1:n, 1, prob = pr)
  }
  kmeans(X, X[C, ])
}

# Here distmat(a, b) should return the distances between the rows of two 
# matrices a and b There may be several implementations in R, one is distmat()
# in package pracma.
#
# Please note that AFAIK it is not clear whether the approach of kmeans++ is 
# really better than, e.g., kmeans with several restarts.

# suppressMessages(require(ggfortify, quietly = T))
# ?autoplot
# suppressMessages(require(useful, quietly = T))
# ?plot.kmeans

suppressMessages(require(factoextra, quietly = T))
library(help = factoextra)

# Plot
?fviz_cluster
  
ddmat <- as.matrix(dd)
km3 <- kmpp(ddmat, k = 3)
fviz_cluster(km3, data = ddmat, 
             geom = "point",
             ellipse.type = "norm",
             palette = "Dark2",
             ggtheme = theme_minimal())

km4 <- kmpp(ddmat, k = 4)
fviz_cluster(km4, data = ddmat, 
             geom = "point",
             ellipse.type = "norm",
             palette = "Dark2",
             ggtheme = theme_minimal())

# Selecting k...
?fviz_nbclust

fviz_nbclust(ddmat, kmeans, method = "wss")

fviz_nbclust(ddmat, kmeans, method = "silhouette")

fviz_nbclust(ddmat, kmeans, method = "gap_stat")

# k = 2, 3, or 4 seem promising.